home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 19 / CU Amiga Magazine's Super CD-ROM 19 (1998)(EMAP Images)(GB)[!][issue 1998-02].iso / CUCD / Programming / LEDA / source / src / arith / igcd.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-16  |  6.8 KB  |  391 lines
  1. /*******************************************************************************
  2. +
  3. +  LEDA  3.1c
  4. +
  5. +
  6. +  igcd.c
  7. +
  8. +
  9. +  Copyright (c) 1994  by  Max-Planck-Institut fuer Informatik
  10. +  Im Stadtwald, 6600 Saarbruecken, FRG     
  11. +  All rights reserved.
  12. + 
  13. *******************************************************************************/
  14.  
  15.  
  16. /*    igcd.c        RD, 21.03.90        */
  17.  
  18. /*    RD, 11.9.91, neues Ibgcd, Verdopplung der Geschwindigkeit. */
  19.  
  20. #include <LEDA/impl/iint.h>
  21. #include <LEDA/impl/iloc.h>
  22.  
  23. void Idgcd(d, a, b)
  24.     Integer *d;
  25.     const Integer *a, *b;
  26. /* d = gcd(a, b);   Euklidischer Algorithmus, Division mit Rest */
  27. {    Integer aa, bb, q, r;
  28.     register pInteger pa, pb, pq, pr, copy;
  29.     if (Ieq0(a)) {
  30.         IasI(d, b);
  31.         d->sign=PLUS;
  32.         return;
  33.     }
  34.     if (Ieq0(b)) {
  35.         IasI(d, a);
  36.         d->sign=PLUS;
  37.         return;
  38.     }
  39.     pa=&aa; pb=&bb; pq=&q; pr=&r;
  40.     cIasI(pa, a);
  41.     cIasI(pb, b);
  42.     cI(pq);
  43.     cI(pr);
  44.     pa->sign=PLUS;
  45.     pb->sign=PLUS;
  46.     while (TRUE) {
  47.         uIdiv(pq, pr, pa, pb);
  48.         if (Ieq0(pr))
  49.             break;
  50.         copy=pa; pa=pb; pb=pr; pr=copy;
  51.     }
  52.     IasI(d, pb);
  53.     dI(pa); dI(pb); dI(pq); dI(pr); 
  54. }        /* Idgcd */
  55.  
  56. #ifdef USE_OLD_BGCD
  57.  
  58. void Ibgcd(d, a, b)
  59.     pInteger d, a, b;
  60. /* d = gcd(a, b);   binaerer Algorithmus */
  61. {    register int toshift=0;
  62.     Integer aa, bb;
  63.     register pInteger pa, pb, swap;
  64.     if (Ieq0(a)) {
  65.         IasI(d, b);
  66.         d->sign=PLUS;
  67.         return;
  68.     }
  69.     if (Ieq0(b)) {
  70.         IasI(d, a);
  71.         d->sign=PLUS;
  72.         return;
  73.     }
  74.     pa=&aa; pb=&bb;
  75.     cIasI(pa, a);
  76.     cIasI(pb, b);
  77.     pa->sign=PLUS;
  78.     pb->sign=PLUS;
  79.     while (Ieven(pa) && Ieven(pb)) {
  80.         toshift++;
  81.         Isr1(pa);
  82.         Isr1(pb);
  83.     }
  84.     while (Ieven(pa)) {
  85.         Isr1(pa);
  86.     }
  87.     while (Ieven(pb)) {
  88.         Isr1(pb);
  89.     }
  90.     while (TRUE) {
  91.         if (IgtI(pb, pa)) {
  92.             swap=pa; pa=pb; pb=swap;
  93.         } else {
  94.             ImiasI(pa, pb);
  95.             if (!Ieq0(pa))
  96.                 while (Ieven(pa))
  97.                     Isr1(pa);
  98.             else
  99.                 break;
  100.     }    }
  101.     IasIslint(d, pb, toshift);
  102.     dI(pa);
  103.     dI(pb);
  104. }        /* Ibgcd */
  105.  
  106. #else
  107.  
  108. void Ibgcd(d, a, b)
  109.     Integer *d;
  110.     const Integer *a, *b;
  111. /*
  112.  * d = gcd(a, b);   binaerer Algorithmus, 
  113.  * direkt auf irv.c aufgesetzt.
  114.  * Achtung:
  115.  * Sind die Zahlen a, b sehr unterschiedlich gross, so ist
  116.  * Ibgcd relativ langsam. Deshalb ist es sinnvoll, vorher
  117.  * eine Division mit Rest durchzuf"uhren.
  118.  *
  119.  */
  120. {    register int l, al, bl;
  121.     int toshift=0;
  122.     Integer aa, bb;
  123.     register PLACE *lp, *avec, *bvec;
  124.     if (!a->length) {
  125.         IasI(d, b);
  126.         d->sign=PLUS;
  127.         return;
  128.     }
  129.     if (!b->length) {
  130.         IasI(d, a);
  131.         d->sign=PLUS;
  132.         return;
  133.     }
  134.     cI(&aa);
  135.     cI(&bb);
  136.     /* Nun zuerst eine Division mit Rest */
  137.     if (IgtI(a, b)) {
  138.         IasIreI(&aa, a, b);
  139.         if (!aa.length) {
  140.             IasI(d, b);
  141.             d->sign=PLUS;
  142.             dI(&aa);
  143.             dI(&bb);
  144.             return;
  145.         }
  146.         IasI(&bb, b);
  147.     } else {
  148.         IasIreI(&bb, b, a);
  149.         if (!bb.length) {
  150.             IasI(d, a);
  151.             d->sign=PLUS;
  152.             dI(&aa);
  153.             dI(&bb);
  154.             return;
  155.         }
  156.         IasI(&aa, a);
  157.     }
  158.     avec=aa.vec;
  159.     bvec=bb.vec;
  160.     al=aa.length;
  161.     bl=bb.length;
  162.     while (!(*avec & 1) && !(*bvec & 1)) {
  163.         toshift++;
  164.         vecsr1(avec, al);
  165.         if (!avec[al-1])
  166.             al--;
  167.         vecsr1(bvec, bl);
  168.         if (!bvec[bl-1])
  169.             bl--;
  170.     }
  171.     while (!(*avec & 1)) {
  172.         vecsr1(avec, al);
  173.         if (!avec[al-1])
  174.             al--;
  175.     }
  176.     while (!(*bvec & 1)) {
  177.         vecsr1(bvec, bl);
  178.         if (!bvec[bl-1])
  179.             bl--;
  180.     }
  181.     while (TRUE) {
  182.         if (bl>al || ((al==bl)&& vecgt(bvec, avec, bl))) {
  183.             lp=avec; avec=bvec; bvec=lp;
  184.             l=al; al=bl; bl=l;
  185.         } else {
  186.             lp=&(avec[al-1]);
  187.                     cvecsubto(avec, bvec, bl);
  188.             while ((al>0)&&(! *lp)) {
  189.                 al--; lp--;
  190.             }
  191.             if (al) 
  192.                 while (!(*avec & 1)) {
  193.                 register PLACE low;
  194.                 l=1;
  195.                 low=*avec>>1;
  196.                 while ( !(low & 1) && l<LOGRADIX-1 ) {
  197.                     l++;
  198.                     low>>=1;
  199.                 }
  200.                 vecsri(avec, al, l);
  201.                 if (!avec[al-1])
  202.                     al--;
  203.                 }
  204.             else
  205.                 break;
  206.     }    }
  207.     if (avec==aa.vec) {
  208.         bb.length=bl;
  209.         bb.sign=PLUS;
  210.         IasIslint(d, &bb, toshift);
  211.         dI(&aa);
  212.         dI(&bb);
  213.         return;
  214.     } else {
  215.         aa.length=bl;
  216.         aa.sign=PLUS;
  217.         IasIslint(d, &aa, toshift);
  218.         dI(&aa);
  219.         dI(&bb);
  220.         return;
  221.     }
  222. }        /* Ibgcd */
  223.  
  224. #endif
  225.  
  226. /*************************************************/
  227.  
  228. void Ielba(d, u, v, a, b)
  229.     Integer *d, *u, *v;
  230.     const Integer *a, *b;
  231. /* d = gcd(a, b) = ua + vb;   Euklid-Lenstra-Berlekamp */
  232. {    Integer aa, bb, wua, wva, wub, wvb, q, r;
  233.     register pInteger pa, pb, pq, pr, ua, va, ub, vb, copy;
  234.     if (Ieq0(a)) {
  235.         IasI(d, b);
  236.         Ias0(u);
  237.         Ias1(v);
  238.         v->sign=b->sign;
  239.         d->sign=PLUS;
  240.         return;
  241.     }
  242.     if (Ieq0(b)) {
  243.         IasI(d, a);
  244.         Ias1(u);
  245.         Ias0(v);
  246.         u->sign=a->sign;
  247.         d->sign=PLUS;
  248.         return;
  249.     }
  250.     pa=&aa; pb=&bb; pq=&q; pr=&r;
  251.     ua=&wua; va=&wva; ub=&wub; vb=&wvb;
  252.     cIasI(pa, a);
  253.     cIasI(pb, b);
  254.     cI(pq);
  255.     cI(pr);
  256.     pa->sign=PLUS;
  257.     pb->sign=PLUS;
  258.     cIasint(ua, 1);
  259.     ua->sign=a->sign;
  260.     cIasint(va, 0);        /* pa == ua*a + va*b */
  261.     cIasint(ub, 0);
  262.     cIasint(vb, 1);
  263.     vb->sign=b->sign;    /* pb == ub*a + vb*b */
  264.     while (TRUE) {
  265.         uIdiv(pq, pr, pa, pb);
  266.         if (Ieq0(pr))
  267.             break;
  268.         copy=pa; pa=pb; pb=pr; pr=copy;
  269.         IasImuI(pr, pq, ub);
  270.         copy=ua; ua=ub; ub=copy;
  271.         ImiasI(ub, pr);
  272.         IasImuI(pr, pq, vb);
  273.         copy=va; va=vb; vb=copy;
  274.         ImiasI(vb, pr);
  275.     }
  276.     IasI(d, pb);
  277.     IasI(u, ub);
  278.     IasI(v, vb);
  279.     dI(pa); dI(pb); dI(pq); dI(pr);
  280.     dI(ua); dI(ub); dI(va); dI(vb); 
  281. }        /* Ielba */
  282.  
  283. void Ibelba(d, u, v, a, b)
  284.     Integer *d, *u, *v;
  285.     const Integer *a, *b;
  286. /* d = gcd(a, b) = u*a + v*b;   binaerer Algorithmus */
  287. {    register int toshift=0;
  288.     Integer aa, bb, A, B, wxa, wxb, wya, wyb;
  289.     register pInteger pa, pb, pA, pB, xa, xb, ya, yb, swap;
  290.     if (Ieq0(a)) {
  291.         IasI(d, b);
  292.         Ias0(u);
  293.         Ias1(v);
  294.         v->sign=b->sign;
  295.         d->sign=PLUS;
  296.         return;
  297.     }
  298.     if (Ieq0(b)) {
  299.         IasI(d, a);
  300.         Ias1(u);
  301.         Ias0(v);
  302.         u->sign=a->sign;
  303.         d->sign=PLUS;
  304.         return;
  305.     }
  306.     pa=&aa; pb=&bb; xa=&wxa; xb=&wxb; ya=&wya; yb=&wyb;
  307.     pA=&A; pB=&B;
  308.     cIasI(pa, a);
  309.     cIasI(pb, b);
  310.     while (Ieven(pa) && Ieven(pb)) {
  311.         toshift++;
  312.         Isr1(pa);
  313.         Isr1(pb);
  314.     }
  315.     cIasI(pA, pa);
  316.     cIasI(pB, pb);
  317.     pa->sign=PLUS;
  318.     pb->sign=PLUS;
  319.     cIasint(xa, 1);
  320.     xa->sign=a->sign;
  321.     cIasint(ya, 0);        /* pa == xa*pA + ya*pB */
  322.     cIasint(xb, 0);
  323.     cIasint(yb, 1);
  324.     yb->sign=b->sign;    /* pb == xb*pA + yb*pB */
  325.     while (Ieven(pa)) {
  326.         Isr1(pa);
  327.         if (Ieven(xa) && Ieven(ya)) {
  328.             Isr1(xa);
  329.             Isr1(ya);
  330.         } else {
  331.             IplasI(xa, pB);
  332.             Isr1(xa);
  333.             ImiasI(ya, pA);
  334.             Isr1(ya);
  335.     }    }
  336.     while (Ieven(pb)) {
  337.         Isr1(pb);
  338.         if (Ieven(xb) && Ieven(yb)) {
  339.             Isr1(xb);
  340.             Isr1(yb);
  341.         } else {
  342.             IplasI(xb, pB);
  343.             Isr1(xb);
  344.             ImiasI(yb, pA);
  345.             Isr1(yb);
  346.     }    }
  347.     while (TRUE) {
  348.         if (IgtI(pb, pa)) {
  349.             swap=pa; pa=pb; pb=swap;
  350.             swap=xa; xa=xb; xb=swap;
  351.             swap=ya; ya=yb; yb=swap;
  352.         } else {
  353.             ImiasI(pa, pb);
  354.             ImiasI(xa, xb);
  355.             ImiasI(ya, yb);
  356.             if (!Ieq0(pa))
  357.                 while (Ieven(pa)) {
  358.                 Isr1(pa);
  359.                 if (Ieven(xa) && Ieven(ya)) {
  360.                     Isr1(xa);
  361.                     Isr1(ya);
  362.                 } else {
  363.                     IplasI(xa, pB);
  364.                     Isr1(xa);
  365.                     ImiasI(ya, pA);
  366.                     Isr1(ya);
  367.                 }    }
  368.             else
  369.                 break;
  370.     }    }
  371.     IasIslint(d, pb, toshift);
  372.     IasI(u, xb); IasI(v, yb);
  373.     dI(pa);    dI(pb);    dI(pA);    dI(pB);
  374.     dI(xa);    dI(xb);    dI(ya);    dI(yb);
  375. }        /* Ibelba */
  376.  
  377. /***********************************************/
  378.  
  379. void Ireduce(a, b)
  380.     pInteger a, b;
  381. /* Kuerze gcd von a und b */
  382. {    Integer d, r;
  383.     cI(&d);
  384.     cI(&r);
  385.     Ibgcd(&d, a, b);
  386.     Idiv(a, &r, a, &d);
  387.     Idiv(b, &r, b, &d);
  388.     dI(&d);
  389.     dI(&r);
  390. }
  391.